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Abstract: The goal of this paper is to relate numerical dissipations that are inherited in high 
order shock-capturing schemes with the onset of wrong propagation speed of discontinuities. For 
pointwise evaluation of the source term, previous studies indicated that the phenomenon of wrong 
propagation speed of discontinuities is connected with the smearing of the discontinuity caused by 
the discretization of the advection term. The smearing introduces a nonequilibrium state into the 
calculation. Thus as soon as a nonequilibrium value is introduced in this manner, the source term 
turns on and immediately restores equilibrium, while at the same time shifting the discontinuity 
to a cell boundary. The present study is to show that the degree of wrong propagation speed 
of discontinuities is highly dependent on the accuracy of the numerical method. The manner 
in which the smearing of discontinuities is contained by the numerical method and the overall 
amount of numerical dissipation being employed play major roles. Moreover, employing finite time 
steps and grid spacings that are below the standard Courant-Friedrich-Levy (CFL) limit on shock- 
capturing methods for compressible Euler and Navier-Stokes equations containing stiff reacting 
source terms and discontinuities reveals surprising counter-intuitive results. Unlike non-reacting 
flows, for stiff reactions with discontinuities, employing a time step and grid spacing that are below 
the CFL limit (based on the homogeneous part or non-reacting part of the governing equations) 
does not guarantee a correct solution of the chosen governing equations. Instead, depending on 
the numerical method, time step and grid spacing, the numerical simulation may lead to (a) the 
correct solution (within the truncation error of the scheme), (b) a divergent solution, (c) a wrong 
propagation speed of discontinuities solution or (d) other spurious solutions that are solutions of the 
discretized counterparts but are not solutions of the governing equations. The present investigation 
for three very different stiff system cases confirms some of the findings of Lafon & Yee (1996) and 
LeVeque & Yee (1990) for a model scalar PDE. The findings might shed some light on the reported 
difficulties in numerical combustion and problems with stiff nonlinear (homogeneous) source terms 
and discontinuities in general. 

Keywords :VLigh order numerical methods, Numerical combustion, Chemical reacting flows, 
Nonequilibrium flows, Numerical methods for stiff source terms with shocks. 


1 Introduction 

Consider 3D reactive Euler equations of the form 

U t + F(U) x + G(U) v + H(U) z = S(U), (1) 
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where U, F(U), G(U), H(U) and S(U) are vectors. Here, S(U) is restricted to be homogeneous in U; that 
is, ( x,y,z ) and t do not appear explicitly in S(U). If the time scale of the ordinary differential equation 
(ODE) Ut = S(U) for the source term is orders of magnitude smaller than the time scale of the homogeneous 
conservation law Ut + F(U) X + G(U) y + H{U) Z = 0, then the problem is said to be stiff due to the source 
terms. In combustion or high speed chemical reacting flows the source term represents the chemical reactions 
which may be much faster than the gas flow. This leads to problems of numerical stiffness due to chemical 
reactions. Insufficient spatial/temporal resolution may cause an incorrect propagation speed of discontinuities 
and nonphysical states for standard dissipative numerical methods that were developed for non-reacting flows. 

This numerical phenomenon was first observed by Colella et al. [1] in 1986 who considered both the 
reactive Euler equations and a simplified system obtained by coupling the inviscid Burgers equation with 
a single convection/reaction equation. LeVeque and Yee [2] showed that a similar spurious propagation 
phenomenon can be observed even with scalar equations, by properly defining a model problem with a stiff 
source term. They introduced and studied the simple one-dimensional scalar conservation law with an added 
nonhomogeneous parameter dependent source term 

u t +u x = S(u), (2) 

S(u) = -fm(u - ^)(u - 1), (3) 

where the parameter ^ can be described as the reaction time. When fj is very large, a wrong propagation 
speed of discontinuity phenomenon by dissipative numerical methods will be observed in coarse grids. In 
order to isolate the problem, LeVeque and Yee solved (2) and (3) by the fractional step method. For this 
particular source term, the reaction (ODE) step of the fractional step method can be solved exactly. In 
their study using pointwise evaluation of the source term ( S(u ) is evaluated at the j grid point index, i.e. , 
S(uj) for each time evolution), the phenomenon of wrong propagation speed of discontinuities is connected 
with the smearing of the discontinuity caused by the spatial discretization of the advection term. They 
found that the propagation error is due to the numerical dissipation contained in the scheme, which smears 
the discontinuity front and activates the source term in a nonphysical manner. The smearing introduces 
a nonequilibrium state into the calculation. Thus as soon as a nonequilibrium value is introduced in this 
manner, the source term turns on and immediately restores equilibrium, while at the same time shifting 
the discontinuity to a cell boundary. By increasing the spatial resolution by an order of magnitude, they 
were able to improve towards the correct propagation speed. It is remarked here that in a general stiff 
source term problem, a sufficient spatial resolution is as important as temporal resolution when the reaction 
step of the fractional step method cannot be solved exactly. As will be shown in the present study, on one 
hand, the degree of wrong propagation speed of discontinuities is highly dependent on the accuracy of the 
numerical method. On the other hand, the manner in which the smearing of discontinuities is contained 
by the numerical method and the overall amount of numerical dissipation being employed play major roles. 
Moreover, employing finite time steps and grid spacings that are below the standard Courant-Friedrich-Levy 
(CFL) limit on shock-capturing methods for compressible Euler and Navier-Stokes equations containing stiff 
reacting source terms and discontinuities reveals surprising counter-intuitive results. 

Based on the work of [2, 3, 4, 5, 6, 7], in addition to the incorrect propagation speed of discontinuities, 
other spurious numerics, that are directly tied to the amount of numerical dissipation contained in the chosen 
scheme and the numerical treatment of source terms may result in 

• Possible spurious steady-state numerical solutions and spurious standing waves [3, 4, 5, 6]: 

It was shown in Lafon & Yee [5, 6] and Griffiths et al. [3] that various ways of discretizing the nonlinear 
reaction terms can affect the stability of, and convergence to, the spurious numerical steady states 
and/or the exact steady states. Pointwise evaluation of the source terms appears to be the least stable. 
The studies of Lafon & Yee [5, 6] indicated that numerical phenomena of incorrect propagation speeds of 
discontinuities may be linked to the existence of some stable spurious steady-state numerical solutions. 
More importantly, the different combination of time step, grid spacing and initial condition plays a 
major role in obtaining the correct solution. In addition, it was shown in Yee et al. [3] and Griffiths 
et al. [4] that spurious discrete traveling waves can exist, depending on the method of discretizing the 
source term. Recently, Wang et al. [8] indicated that a well-balanced scheme for reacting flows can 
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minimize certain spurious numerics. 

Studies linking spurious numerical standing waves for (2) and (3) by first and second-order spatial and 
temporal discretizations can be found in Lafon and Yee [5, 6] and Griffiths, Stuart and Yee [4, 9]. 

• Possible wrong prediction of transition point Reynolds number by DNS due to spurious 
bifurcation that created a false transition point: Inaccuracy of the scheme or insufficient grid 
points might lead to possible spurious bifurcation as well as creating wrong propagation speed of 
discontinuities and smearing of turbulent fluctuations. See [9] for a discussion. 

The term “spurious (numerical) solutions” here refers to computed solutions that are solutions of the dis- 
cretized counterparts but are not solutions of the considered governing equation. Pointwise evaluation of the 
source term here means that, for each time evolution, S(U) is evaluated at the single grid point S(Uj t k,i), 
where (j, k. 1) is the grid point index. 

For the last two decades, the wrong speed phenomenon has attracted a large volume of research work in 
the literature (see, e.g., [10, 11, 12, 13, 14, 9, 15, 16, 17, 18, 19, 20, 21]). Various strategies have been proposed 
to overcome this wrong speed difficulty for one to two species cases with a single reaction. Since numerical 
dissipation that spreads the discontinuity front is the cause of the wrong propagation speed of discontinuities, 
a natural strategy is to avoid any numerical dissipation in the scheme. In combustion, level set and front 
tracking methods were used to track the wave front to minimize this spurious behavior [15, 16, 17, 18]. See 
Wang et al. [22] for a comprehensive overview of the last two decades of development. Wang et al. also 
proposed a new high order finite difference method with subcell resolution for advection equations with stiff 
source terms for a single reaction for (1.1) to overcome the difficulty. Research for multi-species (3 or more 
species and multi-reactions) is forthcoming. 

1.1 Objective and Outline 

This is a follow on work to Wang, Shu, Yee & Sjogreen, Yee, Kotov & Sjogreen and related earlier work [22, 23, 
3, 4, 5, 6[. The objective of this paper is to study spurious behavior of high order shock-capturing methods 
using the pointwise evaluation of stiff homogenous source terms for problems containing discontinuities. 
Pointwise evaluation is used in the current study in spite of the fact that Lafon & Yee [5, 6] and Griffiths 
et al. [4] indicated two decades ago that pointwise evaluation of the source term (for first and second-order 
schemes) appears to be the least stable. They suggested using non-pointwise evaluation of the source term 
that is more compatible with the convection difference operator. The current study presents a more in-depth 
understanding of the pointwise evaluation approach as the majority of the schemes in use for numerical 
combustion and problems containing stiff sources and discontinuities employ this approach. In addition, 
spurious behavior in this type of highly nonlinear coupling system cases using finite time steps and grid 
spacings is not fully understood. 

Special focus is on the behavior of the recently developed finite difference method with subcell resolution 
[22], and the filter counterparts [24, 25] of the high order subcell resolution method as time step and grid 
spacing are refined. The study also accounts for the scheme behavior as the stiffness of the source term 
increases. Early and less extensive study on the subject has been reported in [26]. Comparison with the 
performance of the Harten & Yee second-order TVD method [27, 28], and standard fifth-order and seventh- 
order WENO schemes (WEN05 and WEN07) [29] are included. Although the subcell resolution idea and its 
filter counterparts are applicable to any high order shock-capturing method, here the study is focused on the 
class of WENO schemes. From here on, the subcell resolution counterparts of WEN05 and WEN07 will be 
denoted by WEN05/SR and WEN07/SR, whereas their filter counterparts will be denoted by WEN05fi/SR 
and WEN07fi/SR. 

The outline of this paper is as follows: A practical stiff hypersonic chemical nonequilbrium viscous 
computation is illustrated in Section 2 to motivate the current study. The high order methods with subcell 
resolution and their filter counterparts [22, 24, 25] are summarized in Section 3. The problem setup for the 
two stiff detonation test cases with numerical results comparing the performance among WEN05, WEN07, 
and the associated filter version of WENO 5 (WEN05fi) [24, 25], WEN05/SR and WEN05fi/SR are then 
presented in Section 4. The present investigation for three very different system cases confirms the findings 


3 



of Lafon & Yee and LeVeque & Yee for a model scalar PDE. In all of the computations, the classical fourth- 
order Runge-Kutta method (RK4) and the Roe flux with Roe’s average states [30] are used. Performance 
using the third-order TVD Runge-Kutta [31] is similar but with a slightly smaller CFL limit. All the WENO 
schemes are the original form of Jiang & Shu [29], except for one case where the finite difference form of the 
recently developed positive WENO scheme [32] using the Lax-Friedrichs was tested. 


2 Motivation: An Unsteady Nonequilibrium Navier-Stokes Com- 
putation [26] 

In general, the reacting terms that arise from nonequilibrium flows in hypersonic aeronautics are less stiff 
than their counterparts in combustion. However, there are stiff chemical nonequilibrium flows that are due 
to the reaction terms. Before the study of two stiff detonation test cases, a stiff 13-species, one temperature 
nonequilibrium model related to the NASA Ames Electric Arc Shock Tube (EAST) experiment is briefly 
investigated. Detailed study will be reported in a forthcoming publication [26] . See [33] for a brief introduc- 
tion and earlier simulations. The reason for this introductory example is to illustrate that it is unlike earlier 
work in [2, 3, 4, 5, 6], where detailed analysis using dynamical system theory were possible. A complex high 
Mach number and high temperature problem like EAST is very costly even for a 3D coarse grid complete 
unsteady simulation. The length of the EAST shock tube experiment is very long and the associated flow 
physics is multiscale with multi-reaction terms [33]. 


2.1 Governing Equations 

In component form of (1.1), a 3D nonequilibrium Navier-Stokes system for the 8.5 m (meter) EAST problem 
(with the thermo-nonequilibrium part neglected) for a preliminary study is given by: 


dp s 

dt 


d 

T v: (p s Uj T Psd S j) — 


dxj 


d d 

^{pUi) + Q^ipUiUj +P0ij - T ij) = o 


dt E 


dxj 


Uj(E+p) + qj + E Psdgjha 'll; T,;- 


= 0, 


(4) 

(5) 

(6) 


where U = (p s , pu;, E) are the conservative variables, p s are the partial densities with k = 1,...,N S for a 
mixture of N s species. Here i = 1,2,3 for 3D. u^i = 1,2.3 are the mixture x, y and ^-velocities, E is the 
mixture total energy per unit volume, p is the pressure, K(T) is the chemical reaction rate and T is the 
temperature. The mixture total density, the pressure and the total energy per unit volume are 


N s 


5 — 1 


p = Y J Ps, P = ^ t E^’ E = Y J Ps(e s {T) + h° s 


N s 


-pv 


S— 1 


(7) 


where R is the universal gas constant, hP s are the species formation enthalpies, and M s indicates the species 
molar masses. 

The viscous stress tensor is given by: 


Tij = p 



duj_ 

dxi 


2 du k „ 
^3 dx k 6ij ' 


(8) 


The diffusion flux is given by: 


, n dX s 
SJ ~ s dxi 


where D s is the diffusion coefficient and X s is the mole fraction of species s. 


(9) 
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p 

1.10546 kg/m 6 

T 

6000 K 

P 

12.7116 MPa 

Yn e 

0.9856 

Yn 2 

0.0144 


P 

3.0964 x 10 4 kg/m 6 

T 

300 K 

P 

26.771 Pa 

Yq 2 

0.21 

Yn 2 

0.79 


Table 1: High (left) and low (right) pressure region initial data 


The conduction heat flux is given by: 


Qj 


, = -A 


dT 

dxj' 


where A is the thermal conductivity of the mixture. The chemical source term is given by: 


(10) 


N r 


— M s ^ ) (Ps,r Q's,r') 


N s 

k f,r H 
m = 1 


Pm 

M„ 


N s 

kb,r- JHj 

1 



( 11 ) 


where a and b are the stoichiometric coefficients, and the forward reaction rates kf_ r coefficients are given 
by Arrhenius’ law: 

k Lr = A f)r T n f-r exp {-E f>r /kT). (12) 


The backward reactions rates coefficients are computed as kb, r = kf,r/K^ q r , where K^ q r is the equilibrium 
constant. 


Due to the multiscale and multi-stiffness of the problem [33], numerical simulations for ID (* = 1) and 
2D (* = 1,2) are considered first in [26]. Numerical study of grid size and numerical method dependence 
of the computed shear and shock locations as the grid is refined for ID and 2D simplifications of the 3D 
EAST problem will be illustrated here. All the computations employ a multi- D high order single/overset 
grid nonequilibrium code ADPDIS3D [34]. Due to high computational cost, only non-overset grid results for 
a very early stage of the unsteady flow development are presented. The desired simulation requires that the 
shock wave propagates to a 8.5 meter distance. The MUTATION library [35], developed by Thierry Magin 
and Marco Panesi, is used for the numerical experiment to provide reaction rate and transport properties. 


Here the motivation is to show several standard shock-capturing methods and their filter counterpart 
schemes for the early time evolution of the flow. An earlier study using the subcell resolution idea [22] 
when applied to only one out of the 13 species with at least a dozen stiffness reaction coefficients shows no 
improvement over the standard TVD scheme. See the next subsection for the ID EAST simulation as an 
illustration. Note that the Wang et al. subcell resolution idea was constructed for a single reaction. Here, 
for this viscous simulation, all the CFL values are based on the convection and viscous part of the PDEs. 


2.2 Problem setup for the ID 13 Species EAST Problem 

The computational domain has a total length of 8.5 m. The left part of the domain with length 0.1m is 
a high pressure region. The right part of the domain with length 8.4m is a low pressure region. The gas 
mixture consists of 13 species: 

e“, He, N, O, N 2 , NO, 0 2 , N + , NO + , N+ ,0% , 0 + , He+. 

The initial conditions of the high and low pressure regions are listed in the Table 1. For the left-side boundary 
the Euler (slip) wall condition is applied, and for the right-side, the zero gradient condition is applied for all 
variables. 


2.2.1 13 Species ID EAST Simulation 

Figure 1 shows the computation using the Harten-Yee second-order TVD scheme [27] for three grids: 
501,1001,10001 at time Tend = 0.325 x 10~ 4 s. One can observe the shift in the shear (left discontinu- 
ity) and the shock (right discontinuity) locations as the grid is refined. The width of the space between the 
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Figure 1: 13 species ID EAST problem: Second-order Harten-Yee TVD simulation for three grids: 501, 1001, 10001, 
and Tend = 0.325 x 10~ 4 s, with CFL = 0.8. 


shear and the shock shrinks as the grid is refined. Figure 2 shows the comparison among five methods with 
the reference solution using the TVD scheme with a 10,001 grid. Here 

• TVDafi+split: Sixth-order central base scheme with the Ducros et al. splitting of the governing 
equations. The flow sensor for the filter step is based on the shock and shear locations instead of using 
wavelets. 

• WEN05-llf: WEN05 using the local Lax-Friedrichs flux. 

• WENOSPafi+split: Nonlinear filter counterpart of the positive WEN05 using the local Lax-Friedrichs 
flux. The flow sensor for the filter step is based on the shock and shear locations instead of using 
wavelets. For the finite difference form of the positive WEN05-llf, see [32]. 

• TVD/SR: Finite difference scheme with subcell resolution (on one of the reaction coefficients) using 
the Harten & Yee TVD scheme as the convection difference operator in the fractional step method. 

Figure 2 indicates that the least dissipative scheme (among the considered schemes), the better the 
prediction of the shear and shock locations when compared with the reference solution using the TVD 
scheme with a 10, 001 grid. The result indicates that TVDafi+split is more accurate than WEN05-llf. This 
is due to the fact that TVDafi+split reduces the amount of numerical dissipation away from high gradient 
regions. The TVD/SR using only one stiffness coefficient in the subcell resolution stage gives no improvement 
over the standard TVD (with no subcell resolution). 

2.3 Problem setup for a 13 Species 2D EAST Problem 

The computational domain is half of the 2D shock tube y-height with total length 8.5 m, height 0.0508m and 
symmetry boundary condition imposed on the top. The left part of the a;-domain with length 0.1m is a high 
pressure region. The right part of the domain with length 8.4m is a low pressure region. The gas mixture 
consists of the same 13 species as the ID simulation: 

e~,He, N , O, N 2 , NO, 0 2 , N+,NO + , N+, 0+,0+, He + . 

The initial conditions of the high and low pressure regions are listed in the Table 1. 
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Figure 2: 13 species ID EAST problem: Comparison among 5 methods using a 501 grid with CFL = 0.8, and 
Tend = 0.325 x 10~ 4 s. See text for method notation. 
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Figure 3: Schematic of a 13 species 2D EAST problem. 
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Figure 4: 2D 13 species EAST simulation by TVD for CFL = 0.7 and Tend = 10“ 5 s: Top Row - Three 
2 -direction grid refinement 601 x 121, 1201 x 121 and grid clustering between shear and shock in the re- 
direction of 691 x 121. All y-grid use boundary grid stretching with a minimum of A y = 10 -5 . Bottom Row: 
Two x-direction grid refinement 1201 x 121 and grid clustering between shear and shock in the 2 -direction 
of 691 x 121. All y-grid use boundary grid stretching with a minimum of Ay = 5 x 10” 6 . 


grid Nx 

601 

1201 

1201 

691 

691 

refine 

no 

no 

no 

yes 

yes 

min h y 

l.e-5 

l.e-5 

5.e-6 

l.e-5 

5.e-6 

Shock T max , K 

15,846 

18,851 

18,848 

25,098 

25,015 

Shear T^dcc-) K 

11,301 

11,203 

11,203 

10,598 

10,598 


Table 2: Shock and Shear maximum temperature grid dependence at time Tend = 10 _5 s. N x indicates the 
grid spacing in the ^-direction. The last two columns are for the grid clustering results for two different 
minimum y-grid stretching. 


For the left boundary the slip (Euler) wall condition is applied. For the right-side the zero gradient 
condition is applied for all variables. The bottom boundary is treated as an isothermal wall with the 
constant temperature T. wa u = 300 K. The top boundary is treated as a symmetric boundary condition. 
Figure 3 shows the schematic of the 2D EAST simulation. 

2.3.1 13 Species 2D EAST Simulation 

For this 2D test case a very accurate reference solution is not practical to obtain due to the CPU intensive 
nature of the problem. Here, three levels of refinement are conducted. Figure 3 shows the schematic of 
the 2D EAST simulation at time Tend = 10~ 5 s using CFL = 0.7 by TVD. Figure 4 shows the computed 
temperature contour results by TVD for three levels x- and y-direction grid refinement simulations. The top 
row shows three axdirection grid refinements of 601 x 121, 1201 x 121 and grid clustering between shear and 
shock in the ^-direction of 691 x 121. All y grids use boundary grid stretching with a minimum of Ay = 10 -5 . 
The bottom row shows the same two ^’-direction grid refinements 1201 x 121 and grid clustering between 
shear and shock in the ^-direction of 691 x 121. All y grids use boundary grid stretching with a minimum of 
Ay = 5 x 10 6 . Comparing the two rows of grid refinement study indicates that by refining the axdirection 
grid with the y-direction the same has a big effect on the locations of the shear/shock. This is due to the fact 
that aside from the boundary layer, the shear and shock are nearly one dimensional. However, comparing 
the last two columns of the grid refinement study indicates that by refining the y-direction grid with the 
^-direction the same has no effect on the locations of the shear/shock, but increases the boundary layer 
prediction. As in the ID EAST simulation, the discontinuity locations shift as the ^’-direction grid is refined. 
The width of the distance between the shear and the shock shrank as the grid was refined. The shear and 
shock strength are also different. Table 2 indicates the maximum shear and contact temperature for each 
set of grids. For the minimum grid stretching of Ay = 10~ 5 , the maximum shear temperature is 11,301A^, 
and maximum shock temperature is 15, 846A" for the 601 x 121 grid. However, the shear and shock strength 
are with maximum shear temperature = 11,203A', and maximum shock temperature = 18,851AT for the 
1201 x 121 grid. For the stretched grid the shear and shock strength are with maximum shear temperature 
= 10, 598A', and maximum shock temperature = 25, 098A'. As we decrease the minimum grid stretching 
to Ay = 5 x 10~ 6 , the shear and shock strength are with maximum shear temperature = 11,203A", and 
maximum shock temperature = 18, 848AT for the 1201 x 121 grid. For the stretched grid the shear and shock 
strength are with maximum shear temperature = 10, 598A^, and maximum shock temperature = 25, 015Ah 
Aside from the different shock/shear locations the result indicated in the last column shows the maximum 
temperature at the shock location is higher than the result indicated in the the middle and the first columns. 


These results indicate that the numerical method and grid dependence of the shear and shock locations 
are related to the stiffness of the source terms. Note that for non-reacting flows, numerical method and grid 
dependence of the solution do not affect the location of the discontinuities, but rather affect the degree of the 
smearing of the discontinuities. The implication of the EAST computation exercise is to illustrate the danger 
of practical numerical simulation for problems containing stiff source terms where there is no reliable means 
of assessing the accuracy of the computed result other than by extreme grid refinement, which is beyond the 
capability of the current super computer. Before the detail numerical study for the two stiff detonation test 
cases, the next section gives a brief description of our recently developed high order shock-capturing method 
with specific numerical dissipation controls [22, 24, 25]. 
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3 Overview of Two Recently Developed High Order Shock-Capturing 
Schemes 

Here only the newly developed high order finite difference method with subcell resolution for advection 
equations with stiff source terms ([22]) in 2D is briefly summarized. The key aspects of the filter counterpart 
of the WENO schemes are included at the end of the section. For simplicity of discussion only 2D reactive 
Euler equations are considered. It is noted that the considered schemes are applicable to 3D reactive flows. 
Although the Wang et al. high order scheme with subcell resolution [22] is only developed for a single 
reaction case, the Yee & Sjogreen and Sjogreen & Yee high order nonlinear filter scheme [24, 25, 36, 37, 8] 
is applicable for any number of species and reactions. The high order nonlinear filter scheme with local flow 
sensor is applied to further control the amount of numerical dissipation being used for turbulence with strong 
shocks. 

3.1 2D Reactive Euler Equations 

Consider a 2D inviscid combustion flow containing two species 


{Pl)t + (PlU)x + (PlV)y 

= K{T)P2 

(13) 

{P2)t + (P2U) X + ( P2V) V 

= -K(T)P2 

(14) 

(, pu) t + ( PU 2 + P)x + ( PUV) V 

= 0 

(15) 

( pv) t + (. PUV) X + (pv 2 + p) y 

= 0 

(16) 

Et + (u(E + p)) x + (v(E + p)) y 

= 0 

(17) 


where p\ is the density of burned gas, p 2 is the density of unburned gas, u and v are the mixture x- and 
y- velocities, E is the mixture total energy per unit volume, p is the pressure, K(T ) is the chemical reaction 
rate and T is the temperature. The pressure is given by 

P = (7- !)(£- ^p{u 2 +v 2 ) -qoP 2 ), (18) 

where the temperature T = p/ p and q 0 is the chemical heat released in the reaction. 

The mass fraction of the unburnt gas is 2 = P 2 I p- The mixture density is p = pi + p 2 - 

The reaction rate K(T ) is modeled by an Arrhenius law 

K(T) = K 0 e X p(^^j, (19) 


where Kq is the reaction rate constant and Ti gn 
modeled in the Heaviside form 


K{T) = 


is the ignition temperature. The reaction rate may be also 


f A' 0 T > T. ign 

\ 0 T < T. ign . 


(20) 


3.2 High Order Finite Difference Methods with Subcell Resolution for Advec- 
tion Equations with Stiff Source Terms 

The general fractional step approach based on Strang-splitting [38] for the 2D reactive Euler equations 
written in vector notation 

U t + F(U) X + G(U) y = S{U ) (21) 

is as follows. The numerical solution at time level t n + 1 is approximated by 


Un+i = A f^\ R(At)A(^p\ U n . 


(22) 
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The reaction operator R is over a time step At and the convection operator A is over At/ 2. The two 
half-step reaction operations over adjacent time steps can be combined to save cost. The convection operator 
A is defined to approximate the solution of the homogeneous part of the problem on the time interval, i.e., 


U t + F(U) X + G{U) V =0, t n <t<t 


— L n+ !• 


(23) 


The reaction operator R is defined to approximate the solution on a time step of the reaction problem: 


dU 

dt. 


— S(U), t n < t < t n+ 1- 


(24) 


Here, the convection operator consists of, e.g., WEN05 with Roe flux and RK4 for time discretization. 
If there is no smearing of discontinuities in the convection step, any ODE solver can be used as the reaction 
operator. However, all the standard shock-capturing schemes will produce a few transition points in the shock 
when solving the convection equation. These transition points are usually responsible for causing incorrect 
numerical results in the stiff case. Thus, a direct application of a standard ODE solver at these transition 
points will create incorrect shock speed. To avoid this, here the Harten’s subcell resolution technique [39] in 
the reaction step is employed. The general idea is as follows. If a point is considered a transition point of the 
shock, information from its neighboring points which are deemed not transition points will be used instead. 
In 2D case we apply the subcell resolution procedure dimension by dimension. Here, U T = (pi, p 2 , pu, pv , E) 
and we select the mass fraction 2 as the stiffness indicator. The algorithm proceeds as follows. 

(1) Use a “shock indicator” to identify cells in which discontinuities are believed to be situated. One can 
use any indicator suitable for the particular problem. Here the minmod-based shock indicator in [39, 40] is 
considered. Identify “troubled cell” Rj in both x- and //-directions by applying the shock indicator to, e.g., 


the mass fraction 2 . Define the cell Rj as troubled in the ^-direction if 
with at least one strict inequality, where 


s ij\ > 




and 


> 


■ , 2+l, < 7 I 


sfj = minmod{z i+ i j - 2 ^, 2 ^ - Zi-i,j}. (25) 

Similarly we can define s 1 /^ , the cell Rj as troubled in the //-direction. 

If Rj is only troubled in one direction, we apply the subcell resolution along this direction. If Rj is 
troubled in both directions, we choose the direction which has a larger jump. Namely, if |s?.| > | .s;'U | , subcell 
resolution is applied along the ^-direction, otherwise it is done along the //-direction. In the following steps 
(2)-(3), without loss of generality, we assume the subcell resolution is applied in the x-direction. Assuming 
Rj is troubled in the x-direction, we apply subcell resolution along the x-direction. 

In a troubled cell identified above, we continue to identify its neighboring cells. For example, we can 
define R+ij as troubled if |sf +1 -| > |sf_i -| and |s? +1 -| > |s* +2 -| and similarly define R-ij as troubled if 
[sf_i j\ > |sf_ 2 -| and |s?_ lj | > jsf +1 -|. If the cell R- S ,j and the cell R+ r ,j ( s,r > 0) are the first good cells 
from the left and the right (i.e., R- s + i,j and R+ r - ip are still troubled cells), we compute the fifth-order 
ENO interpolation polynomials pi_ s j(x) and pi +r j(x) for the cells R~ s j and R+ r ,j, respectively. 

(2) Modify the point values Zij, Rj and ptj in the troubled cell Rj by the ENO interpolation polynomials 


Zij =Pi- Stj (xi;z), Tjj = Pi- s ,j(xi;T), = Pi- 8 ,j(xi] p), if 6 > 

%ij = Pi-\-r,j(Xi] z), Tij — (Xii - 0 ? Pij = Pi-\-r,j{x%i P) i ^ 0 X{ 

where the location 6 is determined by the conservation of energy E 


(26) 


/ pi- s j(x;E)dx + / 
'**- 1/2 JO 


i + 1/2 


Pi+ r j(x ; E)dx = EijAx. 


(27) 


Under certain conditions, it can be shown that there is a unique 6 satisfying Eq. (27), which can be solved 
using, for example, a Newton’s method. If there is no solution for 9 or there is more than one solution, we 
choose Rj = Zi+ r j, Rj = R+ r p and pij = pi+ r j- For particular problems one can choose any other suitable 
method for the reconstruction. 
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(3) Use Uij instead of U l;} in the ODE solver if the cell /, ; is a troubled cell. For simplicity, explicit Euler is 
used as the ODE solver. 

(pz)l +1 - ( pz )?• + A (28) 

Here we would like to remark that, implicit temporal discretization cannot be used in this step because the 
troubled values need to be modified explicitly. However, there is no small time step restriction in the explicit 
method used here, because once the stiff points have been modified, the modified source term S(Tjj , Pij , Zij) 
is no longer stiff. Therefore, a regular CFL number is allowed in the explicit method. (Note that if however, 
a linearized form of a two-level implicit time discretization might be suitable for the reaction step operator. 
This will be investigated in the future.) 

Earlier study reported in [22], in general, a regular CFL = 0.1 using the explicit Euler to solve the 
reaction operator step can be used in the subcell resolution scheme to produce a stable solution. But the 
solution is very coarse in the reaction zone because of the underresolved mesh in time. In order to obtain 
more accurate results in the reaction zone, we evolve one reaction step via N r sub steps, i.e., 



in some numerical examples studied in [22] . 

3.3 Well-Balanced High Order Filter Schemes for Reacting Flows ([24, 25, 36, 
37, 8]) 

The high order nonlinear filter scheme of [24, 25, 36] , if used in conjunction with a dissipative portion of 
a well-balanced shock-capturing scheme as the nonlinear numerical flux, is a well-balanced scheme [8[. The 
well-balanced high order nonlinear filter scheme for reacting flows consists of three steps. 

3.3.1 Preprocessing Step 

Before the application of a high order non-dissipative spatial base scheme, the pre-processing step to improve 
stability had split inviscicl flux derivatives of the governing equation(s) in the following three ways, depending 
on the flow types and the desire for rigorous mathematical analysis or physical argument. 

• Entropy splitting of [41] and [42, 43]: The resulting form is non-conservative and the derivation is 
based on entropy norm stability with boundary closure for the initial value boundary problem. 

• The system form of the Ducros et al. splitting [44] : This is a conservative splitting and the derivation 
is based on physical arguments. 

• Tadmor entropy conservation formulation for systems [45]: The derivation is based on mathematical 
analysis. It is a generalization of Tadmor’s entropy formulation to systems and has not been fully 
tested on complex flows. 

See Honein [46] for a comparison of the entropy splitting and other earlier momentum conservation methods. 

3.3.2 Base Scheme Step 

A full time step is advanced using a high order non-dissipative (or very low dissipation) spatially central 
scheme on the split form of the governing partial differential equations (PDEs). Summation-by-parts (SBP) 
boundary operator [47, 48] and matching order conservative high order free stream metric evaluation for 
curvilinear grids [49] are used. High order temporal discretization such as the third-order or fourth-order 
Runge-Kutta (RK3 or RK4) temporal is used. It is remarked that other temporal discretizations can be used 
for the base scheme step. Numerical experiments only focused on RK4 using Roe’s approximate Riemann 
solver. 
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3.3.3 Post-Processing or Nonlinear Filter Step 

Unlike linear spectral filters for spectral methods or linear compact filters for spatially compact (Pade) 
schemes, here the nonlinear dissipative portion of a high order shock-capturing scheme is employed in con- 
junction with a flow sensor to limit the amount of numerical dissipation being used. After the application of 
a non-dissipative high order spatial base scheme on the split form of the governing equation (s), to further im- 
prove nonlinear stability from the non-dissipative spatial base scheme, the post-processing step of [24, 25, 36] 
nonlinearly filtered the solution by a dissipative portion of a high order shock-capturing scheme with a local 
flow sensor. The flow sensor provides locations and amounts of built-in shock-capturing dissipation that can 
be further reduced or eliminated. The idea of these nonlinear filter schemes for turbulence with shocks is 
that, instead of solely relying on very high order high-resolution shock-capturing methods for accuracy, the 
filter schemes [50, 42, 36, 24, 51] take advantage of the effectiveness of the nonlinear dissipation contained in 
good shock-capturing schemes as stabilizing mechanisms (a post-processing step) at locations where needed. 
The nonlinear dissipative portion of a high-resolution shock-capturing scheme can be any shock-capturing 
scheme. For reacting flow, it is best to employ the dissipative portion of a well-balanced shock-capturing 
scheme. By design, the flow sensors, spatial base schemes and nonlinear dissipation models are standalone 
modules. Unlike standard shock-capturing and/or hybrid shock-capturing methods, the nonlinear filter 
method requires one Riemann solve per dimension per time step, independent of time discretizations. The 
nonlinear filter method is more efficient than its shock-capturing method counterparts employing the same 
order of the respective methods. See [25] for the recent improvements of the work [50, 42, 36, 24] that are 
suitable for a wide range of flow speed with minimal tuning of scheme parameters. For all the computations 
shown, the Ducros et al. splitting is employed. This is due to the fact that for the subject test cases we 
need a robust conservative splitting as the preprocessing step. The subcell resolution approach using the 
fractional step procedure can carry over to the aforementioned filter schemes as well. Some attributes of the 
high order filter approach are: 

• Spatial Base Scheme: High order and conservative (no flux limiter or Riemann solver) 

• Physical Viscosity: Contribution of physical viscosity, if it exists, is automatically taken into consider- 
ation by the base scheme in order to minimize the amount of numerical dissipation to be used by the 
filter step 

• Efficiency: One Riemann solve per dimension per time step, independent of time discretizations (less 
CPU time and fewer grid points than their standard shock-capturing scheme counterparts) 

• Accuracy: Containment of numerical dissipation via a local wavelet flow sensor 

• Well-balanced scheme: These nonlinear filter schemes are well-balanced schemes for certain chemical 
reacting flows [8] 

• Stiff Combustion with Discontinuities: For some stiff reacting flow test cases the high order filter 
scheme is able to obtain the correct propagation speed of discontinuities, whereas the standard high 
order shock-capturing (e.g., WENO) schemes cannot (see the result below) 

• Parallel Algorithm: Suitable for most current supercomputer architectures 

The nonlinear filter counterpart of the subcell resolution method employing, e.g., WEN05 or WEN07 
as the dissipative portion of the filter numerical flux (WEN05fi or WEN07fi) can be obtained in a similar 
manner as in Eq. (3.2) with the convection operator replaces by the nonlinear filter scheme and will be 
denoted by WEN05fi/SR or WEN07fi/SR. 

4 Numerical Results 

Here “coarse grids" means standard mesh density requirement for accurate simulation of typical non-reacting 
flows of similar problem setup. The two well known stiff detonation test cases consist of the Arrhenius ID 
Chapman-Jouguet (C-J) detonation wave [20, 21] and a 2D Heaviside detonation wave [19]. These are 
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Figure 5: ID C-J detonation problem, Arrhenius case for the original stiffness Kq at t = 1.8: Pressure and 
density comparison among three standard shock-capturing methods (TVD, WEN05. WEN07) using 50 
uniform grid points with CFL = 0.05. 




Figure 6: ID C-J detonation problem, Arrhenius case for the original stiffness Kq at t = 1.8: Temperature 
and density comparison among standard high order shock-capturing methods and low dissipative methods 
(WEN05.,WEN05/SR, WEN05fi+split and WEN05fi/SR+split) using 50 uniform grid points with CFL = 
0.05. 
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Figure 7: ID C-J detonation problem, Arrhenius case at t = 1.8: Pressure comparison between the original 
stiffness Kq and 4Ao of the source term computed by WEN05 using 50 uniform grid points. All the CFL 
values for the inviscid simulations are based on the convection part of the PDEs. 




Figure 8: C-J detonation problem, Arrhenius case at t = 1.8: Density comparison between WEN05/SR and 
WEN05fi/SR+split for IOORq and 1000 Ay using 50 uniform grid points with CFL = 0.05. 
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the same two test cases considered in [22]. The considered six schemes are WEN05, “WEN05/SR” (the 
newly developed subcell resolution version of WEN05 [22]), “WEN05fi” (the Yee & Sjogreen nonlinear filter 
version of WEN05 using a local flow sensor to further limit the amount of WEN05 numerical dissipation), 
“WEN05fi+split” (the Ducros et al. splitting of the governing equations [44] of WEN05fi [24, 25] ), and 
"WEN05fi/SR+split " (the nonlinear filter version of WEN05/SR with Ducros et al. splitting of the 
governing equations). All of the five methods use the Roe’s average states. For the temporal discretization 
the classical fourth-order Runge-Kutta method (RK4) is used since the TVD RK3 has lower CFL limit than 
RK4 but with a similar behavior as RK4. The results by RK3 are not considered here. Note that all the 
CFL values for the inviscid simulations are based on the convection part of the PDEs. In addition, the 
computed solutions by WEN05 and WEN05/SR presented here could be slightly different from the results 
presented in [22] due to the minor differences in the formulation of the problem (governing equation; e.g., 
different choice of variables) and the use of cell centered and cell-vertex formulation of numerical schemes. 


4.1 ID Chapman- Jouguet (C-J) Detonation Wave (Arrhenius Case) 

The test case is the ID C-J detonation wave (Arrhenius case) [20, 21]. The initial values consist of totally 
burnt gas on the left-hand side and totally unburnt gas on the right-hand side. The density, velocity, and 
pressure of the unburnt gas are given by p u = 1, u u = 0 and p u = 1. 

The initial state of the burnt gas is calculated from C-J condition: 


where 


p 6 = -6 + (6 2 -c) 1 / 2 , 

(30) 

Pu\Pb{l+ 1) -Pu] 
Pb — > 

IPb 

(31) 

— [Pu^il T ( P/PbPb ) /”]/ Pui 

(32) 

Ub = S C J - ( lPb/rho b ) 1/ 2 , 

(33) 

b = -Pu - Puqoii - l), 

(34) 

>1 + 2(7- l)p M p„5o/(7+ !)• 

(35) 


The heat release qo = 25 and the ratio of specific heats is set to 7 = 1.4. The ignition temperature 
Tig n = 25 and Kq = 16,418. The computation domain is [0,30]. Initially, the discontinuity is located at 
x = 10. At time t = 1.8, the detonation wave has moved to x = 22.8. The reference solution is computed by 
the regular WEN05 scheme with 10,000 uniform grid points and CFL=0.05. 


4.1.1 Initial Study of Scheme Behavior [23] 

Figure 5 shows the pressure and density comparison among the standard TVD, WEN05 and WEN07 using 
50 uniform grid points and CFL = 0.05 for the same stiffness K 0 = 16,418 used in [23]. Figure 6 shows 
the pressure and density comparison among the standard WEN05 scheme, WEN05/SR, WEN05fi and 
WEN05fi+split using 50 uniform grid points. For this particular problem and grid size, all standard TVD 
WEN05 and WEN07 exhibit wrong shock speed of propagation with the lower order and more dissipative 
schemes exhibiting the largest error. WEN05fi+split compares well with WEN05/SR for the computed 
pressure solution. WEN05/SR and WEN05fi+split can capture the correct structure using fewer grid 
points than those in [20] and [21]. A careful examination of the 50 coarse grid mass fraction solutions 
indicates that WEN05fi+split is 0.7 grid point ahead of WEN05/SR at the discontinuity location when 
compared to the reference solution. Since WEN05fi+split is less dissipative than WEN05, the restriction of 
the shock-capturing dissipation using the wavelet flow sensor helps to improve the wrong propagation speed 
of discontinuities without the subcell resolution procedure. It is interesting to see that all of the methods 
(except WEN05) produce oscillatory solutions in the vicinity of the reaction front. This behavior prompted 
us to perform a systematic six levels of uniform grid refinements (200, 400, 800, 1600, 3200 and 6400). As the 
number of grid points increases, this oscillatory behavior in the vicinity of the reaction front becomes more 
pronounced. It appears that for this particular test case the new subcell resolution scheme WEN05/SR, 
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Figure 9: ID C-J detonation problem, Arrhenius case using 50 uniform grid points: Density comparison for 
seven CFL numbers by WEN05/SR (left). Number of grid point away from the reference solution (Err) as 
a function of the CFL number (128 CFL values with 6.316455696 x 10 -3 equal increment) for three stiffness 
coefficients (100-ffo, lOOORo, lOOOORo by WEN05SR. A negative "Err" value indicates the number of grid 
points behind the reference shock solution. For certain values of CFL, divergent solutions might occur that 
are outside the plotting area. See e.g., the red and green negative values of Err. All the CFL values for the 
inviscid simulations are based on the convection part of the PDEs. 


WENOfi and WEN05fi+split only perform well with coarse grids. However, for the more dissipative scheme 
WEN05, as we refine the grid, the computed solution gets closer and closer to the reference solution. The 
spurious oscillation might be contributed by the use of the Roe’s average state without any correction for 
reacting flows. See Jenny et al. [52] and related development for details. 

4.1.2 Scheme Behavior with Increase in Stiffness of the Source Terms 

Figure 7 indicates the behavior of WEN05 for two stiffness coefficients of the reaction rate using 50 grid points 
and CFL = 0.05. As the stiffness of the source term increases, the wrong shock location gets further and 
further away from the reference solution. It seems that the reference solution is independent of the stiffness 
coefficient. Figure 8 indicates that as we increase the stiffness coefficient further, WEN05/SR cannot obtain 
the correct shock speed for 1000-Ko stiffness, whereas WEN05fi/SR+split was able to maintain the correct 
shock speed for this grid. For this problem it is indicated in Bao & Jin [19], (Eq. 4.15) that the shock 
speed depends on the initial condition and 7 has a closed form solution. It appears that the shock location 
is independent of the stiffness coefficient for this problem. We use that formula to judge if the reference 
solution is close to the true shock location. For the original A'o case the distance between the reference 
(10, 000 grid) and the exact solution is 5 points which is 0.025 point on the 50 grid point spacing. Due to 
the high cost of obtaining a closer to the exact reference solution, we consider the current reference solution 
as the reference shock location. For the stiffer cases we also use the reference solution for the Kq (although 
the spike at the detonation front is not the same, the shock location should be within one grid point of the 
coarse grid solution). It would be too costly to obtain a better detonation front spike value for the stiffer 
case as all of the coarse grid solutions are far removed from resolving the detonation front. 

4.1.3 Scheme Behavior as a Function of CFL, Grid Refinement and Stiffness of the Source 
Terms 

The result from Figure 7 prompted us to perform a more systematic study on the spurious numerics for the 
test case. Figure 9 shows the effect of the time steps for seven CFL values that are under the CFL limit (left 
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Figure 10: ID C-J detonation problem, Arrhenius case at t = 1.8: Number of grid points away from the 
reference shock solution (Err) as a function of the CFL number (128 discrete CFL values with 6.316455696 x 
10~ 3 equal increment) for three standard shock-capturing methods using 50, 150, 300 uniform grid points 
(across) and for stiffness I\q, 100-Kq, 1000/\q (top to bottom). See Fig. 9 for additional captions. 


18 



Figure 11: ID C-J detonation problem, Arrhenius case at t = 1.8: Number of grid point away from the 
reference shock solution (Err) as a function of the CFL number (128 discrete CFL values with 6.316455696 x 
10~ 3 equal increment) for three low dissipative shock-capturing methods using 50, 150, 300 uniform grid 
points (across) and for stiffness Kq, IOOAo, IOOOATo (top to bottom). See Fig. 9 for additional captions 
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Figure 12: ID C-J detonation problem, Arrhenius case at t = 1.8: Comparison of the same spatial dis- 
cretization with RK4 and RK3 temporal discretization for three low dissipative shock-capturing methods 
using 150, 300 uniform grid points and for stiffness Kq. 


sub-figure), using 50 grid points and WEN05/SR. The right sub-figure shows the error in terms of the number 
of grid points away from the reference shock location (Err) for three stiffness coefficients Ao, lOOA'o and 
lOOOA'o as the function of 128 discrete CFL values. The 128 discrete CFL values are (0.0001 < CFL < 0.803) 
with 6.316455696 x 10~ 3 equal increment. Here, Err is round down to the nearest integer number. Note that 
the CFL limit for WEN05/SR and its filter counterparts are lower than 0.8 due to the explicit Euler reaction 
step. A negative "Err" value indicates the number of grid points behind the reference shock solution. For 
certain values of CFL, divergent solutions might occur that are outside the plotting area. See e.g., the red 
and green negative values of Err. All the CFL values for the inviscid simulations are based on the convection 
part of the PDEs. As the stiffness coefficient increases, it is more and more difficult to obtain the correct 
shock locations by WEN05/SR. 

Figures 10 and 11 illustrate the error for 128 discrete CFL values for the three standard shock-capturing 
schemes (TVD, WEN05 and WEN07) and the three improved high order shock-capturing schemes (WEN05/SR, 
WEN05fi+split and WEN05fi/SR+split). The study is for three uniform grids 50, 150 and 300 (left to right 
columns in the plot) and the three stiffness coefficients Ao, lOOA'o and lOOOA'o (top to bottom in the plot). 
Results indicated that even for CFL = 0.001 using the original Kq stiffness, TVD, and WEN05, are not able 
to obtain the correct shock location using the three considered grid points and the three stiffness coefficients 
as indicated on the Err plot (Fig. 10). For WEN07 for the three grids with the original stiffness fc 0 , the 
correct shock speed can be obtained for most of the CFL values. As a matter of fact, for larger CFL, it 
performed better than WEN05/SR and its filter counterpart. In additions WEN05 produces less “Err” for 
larger CFL. This again indicates that the more accurate scheme results in a better chance of avoiding the 
wrong shock speed spurious numerics. As the stiffness increases, WEN07 no longer produces the correct 
shock speed by the considered three grids. On the contrary, for certain CFL values the improved high order 
shock-capturing methods for reacting flows, e.g., WEN05/SR, WEN05fi+split and WEN05fi/SR+split, are 
able to obtain the correct shock speed. These time steps (CFL values) that can avoid spurious numerics do 
not have to be very small, but they consist of disjoint segments for the time steps that are within the CFL 
limit. It appears that the special dissipation control exhibits more spurious behavior than WEN07 for the 
original Kq case. In addition, WEN05fi/SR+split performs better for the stiffer cases lOOA'o and 1000 A'o 
than the original A' 0 , whereas WENOSfi+split performs better than WEN05 for larger CFL. 

The current study indicated that using the standard CFL condition for the homogeneous part of the 
PDEs (non-reacting part of the governing equations) does not guarantee a correct solution or the correct 
speed of propagation of discontinuities. A stiff ODEs solver with variable time step control in solving the 
reaction part of the operator using the fractional step approach allows the stiffness of the source terms to 
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Figure 13: LeVeque and Yee linear advection and nonlinear stiff source term test case [2]: Number of grid 
points away from the reference shock solution (Err) as a function of the CFL number (128 discrete CFL 
values between (0.001,8) with 6.291338583 x 10 -3 equal increment) by WEN05 and WEN05/SR using 
50, 150,300 uniform grid points (across) and for stiffness K$, 100/io, lOOOA'o (top to bottom). See Fig. 9 for 
additional captions. 
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come into play. However, as indicated in [7, 53, 9, 54], spurious numerics due to the spatial discretization are 
more difficult to avoid because of the nonlinearity of the source terms. The search for further improvement 
of the aforementioned scheme continues. See further discussion on possible improvement on the source term 
treatment numerical strategy in the subsection after next. 

4.1.4 Scheme Behavior by RK3 and by a Single Scalar PDE Case 

All of the results shown are by RK4 temporal discretization. Figure 12 shows that the RK4 and RK3 
exhibit a similar trend but with slight variation in solution behavior for the ID detonation problem. As 
an illustration, the above behavior of the studied schemes also occurs for a simple scalar case (2) and (3) 
studied by LeVeque and Yee [2] in 1990 for second-order schemes. Figure 13 shows a general trend of the 
scheme behavior by WEN05/SR. However, WEN05 behaves differently from the system test case. In this 
case all the nonlinearity and stiffness contained in the governing equation are due to the source term as the 
convection term in the LeVeque & Yee’s scalar model PDE is linear. It appears that the nonlinearity due to 
the convection terms does not alter the general spurious behavior pattern. 

4.2 Are Pointwise Evaluation of the Source Term and Roe’s Average State 
Appropriate? 

On all of the above numerical computations, the pointwise evaluation of the source term was used. However, 
the studies by Lafon & Yee [5, 6] and Griffiths et al. [4] indicated that pointwise evaluation of the source 
term appears to be the least stable. One approach suggested in Lafon & Yee and Griffiths et al. is to use 
non-pointwise evaluation of the source term that is more compatible with the convection difference operator. 
The non-pointwise evaluation of the source term might improve numerical stability and minimize the wrong 
speed of propagation. In addition, there are studies in the literature showing that using the standard Roe’s 
average state for reacting/multi-phrase flows can create spurious oscillations near the discontinuities. See 
for example Jenny, Muller & Thomann [52] and related later articles. Further investigation along these 
directions is planned. The current investigation is to confirm part of the spurious behavior in the studies by 
Lafon & Yee and Griffiths et al. for system cases. 


4.3 2D Detonation Waves 

This example is taken from ([19]). The chemical reaction is modeled by the Heaviside form with the param- 
eters 

7 = 1.4, q 0 = 0.5196 x 10 10 , K 0 = 0.5825 x 10 10 , T ign = 0.1155 x 10 10 

in CGS units. Consider a two-dimensional channel of width 0.005 with solid walls at the upper and lower 
boundaries. The computational domain is [0,0.025] x [0,0.005]. The initial conditions are 


(p,u,v,p, 


z) = 


{pb,u b ,0,p b ,0), if x<£(y), 
(p u ,u u ,0,p u ,l), if x>£(y), 


(36) 


where 


Z(y) = 


0.004 \y - 0.0025] > 0.001, 

0.005 - | y- 0.0025| j y - 0.0025| < 0.001, 


(37) 


and u u = 0, p u = 1.201 x 10 -3 , p u = 8.321 x 10 5 and u b = 8.162 x 10 4 . Values of p b and p b are defined by 
Eq. (30) and (31). In this case u b is greater than defined by Eq. (33). Figure 14 shows the schematic of the 
2D detonation problem. 


4.3.1 Initial Study of Scheme Behavior 

One important feature of this solution is the appearance of triple points, which travel in the transverse 
direction and reflect from the upper and lower walls. A discussion of the mechanisms driving this solution 
is given in [55]. Again, a pointwise evaluation of the source is employed for the 2D test case. Figures 15 and 
16 show the density comparison among the standard WEN05 scheme, WEN05/SR and WEN05fi+split 
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Figure 15: 2D detonation problem at t = 0.3 x 10“ ' and Kq = 0.5825 x 10 10 : Density computed by different 
methods. From left to right: reference solution by the standard WEN05 method using 4000 x 800 uniform 
grid points, WEN05, WEN05/SR and WEN05fi+split using 500x 100 uniform grid points with CFL = 0.05 
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Figure 16: 2D detonation problem at t = 1.7 x 10 - ' and I\q = 0.5825 x 10 10 : Density computed by 
different methods. From left to right: reference solution by the standard WEN05 method using 4000 x 800 
uniform grid points, WEN05, WEN05/SR and WEN05fi+split using 500 x 100 uniform grid points with 
CFL = 0.05. 
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Figure 17: ID cross-section of density at t = 1.7 x 10 -7 by four high order shock-capturing methods for the 
2D detonation problem using 200 x 40 uniform grid points, CFL = 0.05 and Kq = 0.5825 x 10 10 . The left 
figure zoomed in the vicinity of the discontinuity. 


using 500 x 100 uniform grid points at two different times for stiffness A'o = 0.5825 x 10 10 . Figure 17 shows 
the density comparison among the standard WEN05 scheme, WEN05/SR, WEN05fi and WEN05fi+split 
using 200 x 40 and 500 x 100 uniform grid points. The reference solutions are computed by standard WEN05 
with 4000 x 800 grid points. Again, WEN05/SR and WEN05fi+split are able to obtain the correct shock 
speed with similar accuracy. WEN05fi gives a slightly oscillatory solution near x = 0.004. WEN05 and 
WEN05/SR produce no oscillations at the same location. Further improvement of the flow sensor of the 
filter scheme is needed in order to remove the spurious oscillations. Furthermore, for the 500 x 100 grid, 
WEN05fi also obtained the correct shock speed. For CFL = 0.05, however, WEN05fi/SR+split is not able 
to obtain the correct shock speed for the stiff coefficient A'o. 

4.3.2 Scheme Behavior as a Function of CFL, Grid Refinement and Stiffness of the Source 
Terms 

Figure 18 illustrates the error (number of grid points away from the reference shock location) for 128 
discrete CFL values by the three high order shock-capturing schemes WEN05/SR, WEN05fi+split and 
WEN05fi/SR+split. The 128 discrete CFL values are (0.01 < CFL < 0.8) with 6.22047244094488 x 10 -3 
equal increment. For this 2D case, to reduce computational cost, the smallest CFL is 0.01 instead of 0.001 in 
the ID case. Figure 18 shows the error (Err) for two uniform grids 200 x 40 and 500 x 100 (left to right) and 
three stiffness coefficient A' 0 , 100 A'o: 1000A' o (top to bottom). (Note that for the 2D case A' 0 = 0.5825 x 10 10 . 
) Again, as can be seen in this figure, a similar spurious solution behavior as in the ID detonation case 
carries over to the 2D detonation case. However, for this 2D case, WEN05fi+split performs better than 
WEN05fi/SR+split (the reverse of the ID case). Overall, WEN05/SR and WENOSfi+split perform better 
than the other methods. 

4.4 Scheme Performance and Extreme Grid Refinement 

Here, the relative CPU time performance by WEN05/SR, WEN05fi+split and WEN05fi/SR+split using 
the same computer and within the ADPDIS3D code by the pointwise evaluation of the source term is 
included. Fig. 19 shows the ID and 2D detonation problem using 50 uniform grid for CFL = 0.05 and 
RK4 time discretization. In all cases WEN05fi+split and WEN05fi/SR+split consume less CPU time than 
WEN05 and WEN05/SR, respectively. (Note that the larger the number indicated on the table implies 


24 



Figure 18: 2D detonation problem at t = 1.7 x 10~ 7 and K 0 = 0.5825 x 10 10 : Number of grid point 
away from the reference shock solution as a function of the CFL number (128 discrete CFL values with 
6.22047244094488 x 10~ 3 equal increment) for three low dissipative shock-capturing methods using 200 x 40 
and 500 x 100 uniform grid points (across) and for stiffness A'o, 100Ao, lOOOA'o (top to bottom). See Fig. 9 
for additional captions. 
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Scheme Performance 

ID Detonation Problem (Grid 50, CFL = 0.05, RK4) 



WEN05 

WEN05/SR 

WEN05fi+split 

WEN05fi/SR+split 

CPU eff, 
iterations/sec 

1600 

1570 

2050 

1940 

Discontinuity 
location error 
(grid points) 

5 

0 

1 

0 


2D Detonation Problem (Grid 200x40, CFL = 0.05, RK4) 



WEN05 

WEN05/SR 

WEN05fi+split 

WEN05fi/SR+split 

CPU eff, 
iterations/sec 

4.6 

4.3 

10.0 

7.8 

Discontinuity 
location max 

21 

0 

0 

0 

error 

(grid points) 






Figure 19: Sample of scheme performance of WEN05, WEN05/SR, WEN05fi+split and WEN05fi/SR+split for 
CFL = 0.05. 50 grid points are used for the ID case, and 200 x 40 grid points are used for the 2D case with RI<4 as 
the temporal discretization. The CPU times comparison here is based on 8 processor computations. 


less CPU.) Figure 20 shows the extreme refinement computation using 10, 000 grid points for the ID test 
case with CFL = 0.05. It appears that for this particular CFL, WEN05/SR is very close to the reference 
solution but with slight oscillation. WEN05fi/SR+split behaves similarly to WEN05/SR except with an 
increase in small oscillations. However, WEN05fi+split and WEN07fi+ split are not able to obtain the 
correct shock location. This is another counter-intuitive spurious behavior of the considered schemes. 


5 Concluding Remarks 

In [23] we concluded that the filter version of the WEN05 in conjunction with the Ducros et al. splitting 
(WEN05fi+split) is able to obtain the correct propagation speed of discontinuities for two detonation prob- 
lems. The results show that WEN05/SR and WEN05fi+split are able to obtain the correct shock speed 
with similar accuracy, whereas this is not the case for WEN05 & WENOSfi using the same coarse grids. 
Using its original form [25] without further modification, the accuracy of WEN05fi+split was found to be 
nearly as good as WEN05/SR. That conclusion was for one single CFL = 0.05 and the original K 0 stiffness. 

With the present more extensive study the above summary of the scheme behavior needs to be quantified. 
The behavior of these high order shock-capturing schemes is more complicated and does not fall in the 
standard non-reacting flow numerical solution behavior and practices. Aside from the accuracy of the 
scheme, the manner in which the spreading of discontinuities is contained plays a major role in obtaining the 
correct shock location. Choosing the right combination of time step and grid spacing also plays an equal role. 
Several counter-intuitive spurious behaviors are observed as discussed in the numerical result sections. For 
certain instances, smaller CFLs (not extremely small but practical for computation) exhibit more spurious 
behavior. Traditionally, for non-separable finite difference methods, a bigger CFL would give more accurate 
solutions for non-reacting problems, e.g., the MacCormack method. For problems with nonlinear stiff source 
terms, in some instances, it is the larger CFL (within the limit) which exhibits less spurious behavior. The 
results imply that the traditional concept of CFL guideline needs to be revised. Unlike the von Neumann 
analysis for constant coefficient model PDEs containing zero source terms, the linearized stability region 
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Figure 20: ID C-J detonation problem, Arrhenius case at t = 1.8: Behavior of WEN05/SR, WEN05fi+split and 
WEN05fi/SR+split under extreme grid refinement with CFL = 0.05 and 10,000 grid points. The value "k" is the 
k value to control the amount of numerical dissipation indicated in the formula for the filter numerical fluxes [25]. 


for nonhomogeneous PDEs can consist of disjoint intervals, instead of a single continuous interval. The 
implication is that in practical computations where the exact values of these intervals are not known, one 
can easily land in regions that exhibit spurious solutions. One might suspect that our CFL guideline of using 
the homogeneous part of the governing equation is to blame. However, for very small CFL, the stiffness due 
to the reaction term has been accounted for. 

In spite of the counter intuitive results, overall, the more accurate the numerical method, especially the 
less dissipative scheme in conjunction with the containment of spreading the discontinuity, the better the 
performance for very coarse grids (based on fixed grid spacing studies). It performs better than most of the 
previously suggested improved methods reported in the literature for problems containing stiff source terms 
and discontinuities. The subcell resolution method and its nonlinear filter counterparts delay the onset of 
wrong speed of propagation for stiffer coefficients on the same two stiff detonation test cases more than the 
methods reported in the literature. This study also indicated that since this type of scheme is designed for 
coarse grids and moderate stiff source terms, it has additional spurious behavior as the grid is refined and 
the stiffness is further increased. This finding might shed some light on the reported difficulties in numerical 
combustion and problems with stiff nonlinear (homogeneous) source terms and discontinuities in general. 

In order to get a first hand examination of the behavior for practical problems, simplified EAST exper- 
iment setup simulations for a 13 species nonequilibrium flow were conducted. Due to the CPU intensive 
nature of the flow, less in-depth numerical investigations than for the two detonation test cases were con- 
ducted. Results indicate that the numerical method and grid dependence of the shear and shock locations 
are related to the stiffness of the source terms. The reason is that for non-reacting flows, numerical method 
and grid dependent solutions do not affect the location of the discontinuities, but rather change the degree 
of the smearing of the discontinuities. The implication of this exercise is to illustrate the danger of practical 
numerical simulation for problems containing stiff source terms where there is no reliable means of assessing 
the accuracy of the computed result other than by extreme grid refinement as good and reliable experimental 
data are not available . This extreme grid refinement approach is beyond the capability of the current super 
computer for most practical simulations. 

Several thoughts on the causes of the observed spurious behavior that are topics of future research are: (a) 
the spurious oscillations in the vicinity of discontinuities might be due to the use of Roe’s average states [52] , 
(b) the use of a stiff ODE solver with adaptive error control might alleviate some of the spurious numerics 
due to the reaction operator (however, it might present complications in the subcell resolution approach), 
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and (c) as discussed in the ID test case section, the non-pointwise evaluation of the source term that is more 
compatible with the convection difference operator might play a major role in minimizing spurious numerics. 
Studies by Lafon & Yee [5, 6] and Griffiths et al. [4] indicated that pointwise evaluation of the source term 
appears to be the least stable for higher than first-order numerical methods. All three of the above will be 
subjects of the future investigation with emphasis on (c) for higher than second-order methods. 
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